Explore TMB and number of mutations across multiple timepoints of the PBTA Cohort

Tumor evolution project

Data used

In this notebook, we are using the tmb_genomic.tsv file generated from the 01-preprocess-data.Rmd script.

Set up

suppressPackageStartupMessages({
  library(tidyverse)
})

Directories and File Inputs/Outputs

# Detect the ".git" folder. This will be in the project root directory.
# Use this as the root directory to ensure proper sourcing of functions
# no matter where this is called from.
root_dir <- rprojroot::find_root(rprojroot::has_dir(".git"))
scratch_dir <- file.path(root_dir, "scratch")
analysis_dir <- file.path(root_dir, "analyses", "tmb-vaf-longitudinal") 
input_dir <- file.path(analysis_dir, "input")

# Input files
tmb_genomic_file <- file.path(scratch_dir, "tmb_genomic.tsv")
palette_file <- file.path(root_dir, "figures", "palettes", "tumor_descriptor_color_palette.tsv")

# File path to plots directory
plots_dir <-
  file.path(analysis_dir, "plots")
if (!dir.exists(plots_dir)) {
  dir.create(plots_dir)
}

source(paste0(analysis_dir, "/util/function-create-barplot.R"))
source(paste0(analysis_dir, "/util/function-create-dumbbell-plot.R"))
source(paste0(root_dir, "/figures/scripts/theme.R"))

Read in data and process

# Read and process tmb_genomic file
df_total <- readr::read_tsv(tmb_genomic_file, guess_max = 100000, show_col_types = FALSE) %>% 
  group_by(Kids_First_Participant_ID) %>% 
  mutate(cg_distinct = n_distinct(cancer_group) > 1) # to identify samples with different diagnosis across timepoints

# Are there any samples with both WGS and WXS? 
df_total %>% 
  unique() %>% 
  arrange(Kids_First_Participant_ID, experimental_strategy) %>%
  group_by(Kids_First_Participant_ID) %>%
  dplyr::summarise(experimental_strategy_sum = str_c(experimental_strategy, collapse = ";")) 

# There are, so let's remove these from downstream analyses.
df <- df_total %>% 
  filter(!experimental_strategy == "WXS") %>% 
  dplyr::mutate(patient_id = paste(short_histology, Kids_First_Participant_ID, sep = "_")) %>% 
  distinct(cancer_group, .keep_all = TRUE) %>% 
  summarise(cg_sum = str_c(cancer_group, collapse = ",")) %>% # to identify cases with multiple diagnosis
  left_join(df_total, by = c("Kids_First_Participant_ID")) %>% 
  select(Kids_First_Participant_ID, Kids_First_Biospecimen_ID, cg_sum, cancer_group, short_histology, tumor_descriptor, descriptors, tumor_descriptor_sum, tmb, mutation_count)

# How many bs_samples per cg_sum?
print(df %>% count(cg_sum) %>% arrange(desc(n))) 

# A tibble: 37 × 2
   cg_sum                                 n
   <chr>                              <int>
 1 Diffuse midline glioma                52
 2 Ependymoma                            45
 3 Low-grade glioma                      34
 4 High-grade glioma                     28
 5 Medulloblastoma                       26
 6 Atypical Teratoid Rhabdoid Tumor      25
 7 Chordoma                              12
 8 Meningioma                            10
 9 Ewing sarcoma                          7
10 Adamantinomatous Craniopharyngioma     6
# ℹ 27 more rows

# Let's summarize cancer groups with < 10 bs_samples as Other and use this for visualization purposes
cg_sum_df <- df %>% 
  count(cg_sum) %>% 
  dplyr::mutate(cg_sum_n = glue::glue("{cg_sum} (N={n})"))

df <- df %>% 
  left_join(cg_sum_df, by = c("cg_sum")) %>% 
  mutate(cg_plot = case_when(n < 10 ~ "Other",
                             TRUE ~ cg_sum))

# How many bs_samples per cg_plot?
print(df %>% count(cg_plot) %>% arrange(desc(n))) 

# A tibble: 9 × 2
  cg_plot                              n
  <chr>                            <int>
1 Other                               89
2 Diffuse midline glioma              52
3 Ependymoma                          45
4 Low-grade glioma                    34
5 High-grade glioma                   28
6 Medulloblastoma                     26
7 Atypical Teratoid Rhabdoid Tumor    25
8 Chordoma                            12
9 Meningioma                          10

# Read color palette
palette_df <- readr::read_tsv(palette_file, guess_max = 100000, show_col_types = FALSE)

# length(unique(df$Kids_First_Participant_ID))

TMB per Patient case

We will explore TMB per Kids_First_Participant_ID over time by creating stacked barplots.

# Define parameters for function
ylim = max(df$tmb)

# Re-order df
f <- c("Second Malignancy", "Unavailable", "Deceased", "Recurrence", "Progressive", "Diagnosis") # Level df by timepoints
df_plot <- df %>% 
  dplyr:::mutate(tumor_descriptor = factor(tumor_descriptor),
                 tumor_descriptor = fct_relevel(tumor_descriptor, f)) 

# Run function
fname <- paste0(plots_dir, "/", "TMB-genomic-total.pdf")
print(fname)

[1] "/home/rstudio/pbta-tumor-evolution/analyses/tmb-vaf-longitudinal/plots/TMB-genomic-total.pdf"

p <- create_stacked_barplot(tmb_df = df_plot, ylim = ylim)

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

Warning: Removed 1 rows containing missing values (`geom_col()`).

Warning: Removed 1 rows containing missing values (`geom_bar()`).

pdf(file = fname, width = 15, height = 6)
print(p)

Warning: Removed 1 rows containing missing values (`geom_col()`).
Removed 1 rows containing missing values (`geom_bar()`).

dev.off()

png 
  2 

Attention: Hypermutant TMB defined as ≥10 Mb, and Ultrahypermutant TMB defined as ≥100 mutations/Mb in pediatric brain tumors (https://pubmed.ncbi.nlm.nih.gov/29056344/).

Here, we notice that there are samples with high TMB (hyper-mutant samples). Next, we will exclude these samples (threshold >= 10) from downstream analysis. Attention is needed in cases with high number of mutations in only one timepoint as this will lead to un-matched longitudinal samples. We will also remove those so we always have matched longitudinal samples.

# Filter df and remove any samples with single timepoints
df_plot_filter <- df %>%
  filter(!tmb >= 10) %>%
  unique() %>% 
  arrange(Kids_First_Participant_ID, tumor_descriptor) %>%
  group_by(Kids_First_Participant_ID) %>%
  dplyr::summarise(tumor_descriptor_sum = str_c(tumor_descriptor, collapse = ";")) %>% 
  filter(!tumor_descriptor_sum %in% c("Diagnosis", "Progressive", "Recurrence", "Second Malignancy", "Unavailable", "Deceased", "Progressive;Progressive")) %>% 
  dplyr::left_join(df, by = c("Kids_First_Participant_ID", "tumor_descriptor_sum")) %>% 
  mutate(tumor_descriptor = factor(tumor_descriptor),
         tumor_descriptor = fct_relevel(tumor_descriptor, f)) %>% 
  drop_na(tmb) 

# length(unique(df_plot_filter$Kids_First_Participant_ID))

# Define parameters for function
ylim = max(df_plot_filter$tmb)
df_plot_filter <- df_plot_filter

# Run function
fname <- paste0(plots_dir, "/", "TMB-genomic-no-hypermutants.pdf")
print(fname)

[1] "/home/rstudio/pbta-tumor-evolution/analyses/tmb-vaf-longitudinal/plots/TMB-genomic-no-hypermutants.pdf"

p <- create_stacked_barplot(tmb_df = df_plot_filter, ylim = ylim)

Warning: Removed 1 rows containing missing values (`geom_col()`).

Warning: Removed 1 rows containing missing values (`geom_bar()`).

pdf(file = fname, width = 25, height = 8)
print(p)

Warning: Removed 1 rows containing missing values (`geom_col()`).
Removed 1 rows containing missing values (`geom_bar()`).

dev.off()

png 
  2 

TMB across timepoints and cancer types per Patient case

We will explore TMB per cancer group over time by creating dumbbell plots. We classified by using cancer types with the highest number of samples (High- and Low-grade gliomas) versus any other cancer groups.

# How many bs_samples per kids_id and cancer group?
# print(table(df_plot_filter$cg_plot))
print(df_plot_filter %>% 
        count(cg_plot, Kids_First_Participant_ID))

# A tibble: 51 × 3
   cg_plot                          Kids_First_Participant_ID     n
   <chr>                            <chr>                     <int>
 1 Atypical Teratoid Rhabdoid Tumor PT_DVXE38EX                   4
 2 Atypical Teratoid Rhabdoid Tumor PT_HE8FBFNA                   2
 3 Atypical Teratoid Rhabdoid Tumor PT_VTG1S395                   2
 4 Diffuse midline glioma           PT_5NS35B66                   2
 5 Diffuse midline glioma           PT_JNEV57VK                   2
 6 Diffuse midline glioma           PT_JSFBMK5V                   2
 7 Diffuse midline glioma           PT_NK8A49X5                   5
 8 Diffuse midline glioma           PT_P571HTNK                   2
 9 Diffuse midline glioma           PT_PR4YBBH3                   2
10 Ependymoma                       PT_3R0P995B                   2
# ℹ 41 more rows

# Dumbbell plot per cancer group
cancer_groups <- unique(as.character(df_plot_filter$cg_plot))
cancer_groups <- sort(cancer_groups, decreasing = FALSE)
print(cancer_groups)

[1] "Atypical Teratoid Rhabdoid Tumor" "Diffuse midline glioma"          
[3] "Ependymoma"                       "High-grade glioma"               
[5] "Low-grade glioma"                 "Medulloblastoma"                 
[7] "Meningioma"                       "Other"                           

for (i in seq_along(cancer_groups)) {
  print(i)
  df_ct_sub <- df_plot_filter %>% 
    filter(cg_plot == cancer_groups [i])
  
      if (i %in% c(3, 7, 8)) {
    print(cancer_groups [i])
    # Define parameters for function
    ylim <- 2
    } else if (i == 2) {
      print(cancer_groups [i])
      # Define parameters for function
      ylim <- 6
      } else {
        print(cancer_groups [i])
        # Define parameters for function
        ylim <- 4
      }
    

    # Name plots
    fname <- paste0(plots_dir, "/", cancer_groups[i], "-TMB-dumbbell", ".pdf")
    print(fname)
    
    # Run function
    p <- create_dumbbell_ct(tmb_df = df_ct_sub, 
                                 ylim = ylim, 
                                 ct_id = cancer_groups[i])
    pdf(file = fname, width = 12, height = 8)
    print(p)
    dev.off()
}

[1] 1
[1] "Atypical Teratoid Rhabdoid Tumor"
[1] "/home/rstudio/pbta-tumor-evolution/analyses/tmb-vaf-longitudinal/plots/Atypical Teratoid Rhabdoid Tumor-TMB-dumbbell.pdf"

[1] 2
[1] "Diffuse midline glioma"
[1] "/home/rstudio/pbta-tumor-evolution/analyses/tmb-vaf-longitudinal/plots/Diffuse midline glioma-TMB-dumbbell.pdf"

Warning: Removed 1 row containing missing values (`geom_line()`).

Warning: Removed 1 rows containing missing values (`geom_point()`).

Warning: Removed 1 row containing missing values (`geom_line()`).

Warning: Removed 1 rows containing missing values (`geom_point()`).

[1] 3
[1] "Ependymoma"
[1] "/home/rstudio/pbta-tumor-evolution/analyses/tmb-vaf-longitudinal/plots/Ependymoma-TMB-dumbbell.pdf"

[1] 4
[1] "High-grade glioma"
[1] "/home/rstudio/pbta-tumor-evolution/analyses/tmb-vaf-longitudinal/plots/High-grade glioma-TMB-dumbbell.pdf"

[1] 5
[1] "Low-grade glioma"
[1] "/home/rstudio/pbta-tumor-evolution/analyses/tmb-vaf-longitudinal/plots/Low-grade glioma-TMB-dumbbell.pdf"

[1] 6
[1] "Medulloblastoma"
[1] "/home/rstudio/pbta-tumor-evolution/analyses/tmb-vaf-longitudinal/plots/Medulloblastoma-TMB-dumbbell.pdf"

[1] 7
[1] "Meningioma"
[1] "/home/rstudio/pbta-tumor-evolution/analyses/tmb-vaf-longitudinal/plots/Meningioma-TMB-dumbbell.pdf"

[1] 8
[1] "Other"
[1] "/home/rstudio/pbta-tumor-evolution/analyses/tmb-vaf-longitudinal/plots/Other-TMB-dumbbell.pdf"

Total number of mutations across timepoints and biospecimen sample per Patient case

Here, we want to explore the number of mutations per timepoint and biospecimen sample per patient case.

samples <- unique(as.character(df_plot_filter$Kids_First_Participant_ID))
print(samples)

 [1] "PT_02J5CWN5" "PT_0DWRY9ZX" "PT_1ZAWNGWT" "PT_2FVTD0WR" "PT_2MZPGZN1"
 [6] "PT_2YT37G8P" "PT_37B5JRP1" "PT_39H4JN6H" "PT_3GYW6P6P" "PT_3P3HARZ2"
[11] "PT_3R0P995B" "PT_3T3VGWC6" "PT_3VCS1PPF" "PT_5NS35B66" "PT_5ZPPR06P"
[16] "PT_62G82T6Q" "PT_82MX6J77" "PT_89XRZBSG" "PT_962TCBVR" "PT_98QMQZY7"
[21] "PT_99S5BPE3" "PT_B5DQ8FF0" "PT_BRVGRXQY" "PT_BZCJMEX8" "PT_CWXSP19D"
[26] "PT_CXT81GRM" "PT_DFQAH7RS" "PT_DNAJYFZT" "PT_DVXE38EX" "PT_FN4GEEFR"
[31] "PT_HE8FBFNA" "PT_HHG37M6W" "PT_JNEV57VK" "PT_JSFBMK5V" "PT_KMHGNCNR"
[36] "PT_NK8A49X5" "PT_P571HTNK" "PT_PF04R0BH" "PT_PR4YBBH3" "PT_QH9H491G"
[41] "PT_S2SQJVGK" "PT_T4VN7ZRB" "PT_TKWTTRQ7" "PT_TP6GS00H" "PT_TRZ1N1HQ"
[46] "PT_VTG1S395" "PT_XA98HG1C" "PT_XHYBZKCX" "PT_YK7AD0KK" "PT_Z4GS3ZQQ"
[51] "PT_ZMKMKCFQ"

for (i in seq_along(samples)) {
  print(i)
  tmb_sub <- df_plot_filter %>%
    filter(Kids_First_Participant_ID == samples[i])
  
  # Define parameters for function
  ylim = max(df_plot_filter$mutation_count)
 
  # Run function
  p <- create_barplot_sample(tmb_df = tmb_sub,
                             ylim = ylim,
                             sid = samples[i])
  print(p)
}

[1] 1

[1] 2

[1] 3

[1] 4

[1] 5

[1] 6

[1] 7

[1] 8

[1] 9

[1] 10

[1] 11

[1] 12

[1] 13

[1] 14

[1] 15

[1] 16

[1] 17

[1] 18

[1] 19

[1] 20

[1] 21

[1] 22

[1] 23

[1] 24

[1] 25

[1] 26

[1] 27

[1] 28

[1] 29

[1] 30

[1] 31

[1] 32

[1] 33

[1] 34

[1] 35

[1] 36

[1] 37

[1] 38

[1] 39

[1] 40

[1] 41

[1] 42

[1] 43

[1] 44

[1] 45

[1] 46

[1] 47

[1] 48

[1] 49

[1] 50

[1] 51

sessionInfo()

R version 4.2.3 (2023-03-15)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ggthemes_4.2.4  lubridate_1.9.2 forcats_1.0.0   stringr_1.5.0  
 [5] dplyr_1.1.1     purrr_1.0.1     readr_2.1.4     tidyr_1.3.0    
 [9] tibble_3.2.1    ggplot2_3.4.0   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] highr_0.10       bslib_0.4.2      compiler_4.2.3   pillar_1.9.0    
 [5] jquerylib_0.1.4  tools_4.2.3      bit_4.0.5        digest_0.6.31   
 [9] timechange_0.2.0 jsonlite_1.8.4   evaluate_0.20    lifecycle_1.0.3 
[13] gtable_0.3.3     pkgconfig_2.0.3  rlang_1.1.0      cli_3.6.1       
[17] parallel_4.2.3   yaml_2.3.7       xfun_0.38        fastmap_1.1.1   
[21] withr_2.5.0      knitr_1.42       generics_0.1.3   vctrs_0.6.2     
[25] sass_0.4.5       hms_1.1.3        bit64_4.0.5      rprojroot_2.0.3 
[29] tidyselect_1.2.0 glue_1.6.2       R6_2.5.1         fansi_1.0.4     
[33] vroom_1.6.1      rmarkdown_2.21   farver_2.1.1     tzdb_0.3.0      
[37] magrittr_2.0.3   scales_1.2.1     htmltools_0.5.5  colorspace_2.1-0
[41] labeling_0.4.2   utf8_1.2.3       stringi_1.7.12   munsell_0.5.0   
[45] cachem_1.0.7     crayon_1.5.2